function Step5generate_costShocks_17X31()
% *************************************************************************
% GENERATING COST SHOCKS - 17 X 31 (in baseline dimensions: 17 (manufacturing) sectors 31 countries)
% *************************************************************************
%       (1) This code collects the following matrices:
%           - PPI
%           - WIOD related matrices: GAMMA, D, ALPHA, M(SIGMA), SECTORWEIGHTS, TOWEIGHTS (total output weights)
%           - Exchange rate
%       (2) Checks country, sector or country-sector alignments across matrices
%       (3) Generates GAMMA matrices for the following counterfactuals:
%           - autarky, - smooth world I, - smooth world II
%       (4) Generates cost shocks under
%           Baseline: beta_cost = beta_fx = 1, beta_comp = 0
%           Price complementarities: beta_cost = beta_fx = beta, beta_comp = 1 - beta
%           Price complementarities: sector-specific betas: beta_cost = beta_fx = beta_using_sector, beta_comp = 1 - beta_using_sector
%           No exchange rate: beta_cost = beta, beta_fx = beta_comp = 0
%           No FX: E_0 = 0, beta_cost = beta_fx = 1, beta_comp = 0
%           Mechanical pass-through: beta_cost = beta_fx = beta, beta_comp = 0
%           Mechanical pass-through: sector-specific betas: beta_cost = beta_fx = beta_using_sector, beta_comp = 0
%           Input linkages: Balanced I & II, domestic linkages: beta_cost = beta_fx = 1, beta_comp = 0
%       (5) Output is converted to a table & saved in
%           - country-sector level : 'costShocks_allPeriods_table.xlsx'
%           - aggregated to country level : 'aggregatedCostShocks_allPeriods_table.xlsx'
% *************************************************************************
clear; clc;
tic;

%% Settings
rootfolder      = input('Enter the location for the rootfolder in single quotes (''''): \n');
projectfolder   = [rootfolder '\analysis'];
ppifolder       = [projectfolder '\PPI\MatlabFiles\'];
gammafolder     = [projectfolder '\WIOD\GAMMA\MatlabFiles\'];
alphafolder     = [projectfolder '\WIOD\ALPHA\MatlabFiles\'];
sigmafolder     = [projectfolder '\WIOD\SIGMA\MatlabFiles\'];
sweightsfolder  = [projectfolder '\WIOD\SWEIGHTS\MatlabFiles\'];
toweightsfolder = [projectfolder '\WIOD\TOWEIGHTS\MatlabFiles\'];
fxfolder        = [projectfolder '\FX\MatlabFiles\'];
outputfolder    = [projectfolder '\COSTSHOCKS\17X31\'];
betafolder      = [rootfolder '\data\6_BETA\'];
if exist(outputfolder,'dir') == 0
    mkdir(outputfolder);
end
cd(projectfolder);

% Other input files
beta_file = [betafolder 'sectorspecificbetas.xlsx'];
%% -- (Start looping over the year & month)
fprintf([' * Generating cost shocks (17 X 31) for the year: \n']);

for y = 1995:2011
    yr = num2str(y);
    fprintf([' - ' yr ' \n']);

    for m = 1:12
        mo = num2str(m);
        %fprintf([' ************************************** \n  * Running year: ' yr ', month: ' mo '  \n ************************************** \n']);
        
        %% Download vectors & matrices
        
        % Download Gamma matrix
        GAMMAdata = importdata([gammafolder,'WIOD_gamma_' yr '.csv']);
        Gamma_text = GAMMAdata.textdata(2:end,1:3);
        GAMMA = GAMMAdata.data;
        
        % List of countries & sectors
        list_country = unique(Gamma_text(:,2));
        n_country = size(list_country,1);
        list_sector = unique(Gamma_text(:,3));
        n_sector = size(list_sector,1);
        n_cs = n_country * n_sector;
        
        % Generate identity & 1 & 0 matrices
        I = eye(size(GAMMA));
        ONE = ones(size(GAMMA));
        ZERO = zeros(size(GAMMA));
        
        % Generate D matrix
        D = diag(1 - sum(GAMMA,2));
        D_inv = D \ I;
        
        % Download PPI vector
        ppidata = importdata([ppifolder,'PPI_' yr mo '.csv']);
        ppi_text = ppidata.textdata(2:end,1:5);
        PPI = ppidata.data;
        % -size check
        assert(size(PPI,1) == n_cs )
        assert(size(PPI,2) == 1 )
        
        % Download PPI vector
        ppiusddata = importdata([ppifolder,'PPI_inUSD_' yr mo '.csv']);
        ppiusd_text = ppiusddata.textdata(2:end,1:5);
        PPIusd = ppiusddata.data;
        % -size check
        assert(size(PPIusd,1) == n_cs )
        assert(size(PPIusd,2) == 1 )
        
        % Download E matrix - short version
        Edata = importdata([fxfolder,'FX_' yr mo '.csv']);
        E_text = Edata.textdata(2:end,1:4);
        E_short = Edata.data;
        % -size check
        assert(size(E_short,1) == n_country )
        assert(size(E_short,2) == 1 )
        
        % Compute E matrix - Kronecker product: E*1's
        E_0 = kron(E_short,ones(n_sector,1));
        
        % Download Alpha matrix
        ALPHAdata = importdata([alphafolder,'WIOD_alpha_' yr '.csv']);
        Alpha_text = ALPHAdata.textdata(2:end,1:3);
        ALPHA = ALPHAdata.data;
        % -row sum check
        try
            assert(round(sum(sum(ALPHA,1))) == size(ALPHA,2));
        catch
            % Exception: SWE sector 19 row sum equals 0 from 2009
            sumZeroAlpha = Alpha_text(sum(ALPHA,1) == 0,1:3);
            exception = [yr {'SWE'} {'19'}];
            assert(sum(strcmp(sumZeroAlpha,exception)) == 3);
            
            ALPHA(479,479) = 1; % put weight 1 to SWE 19 (itself)
            assert(round(sum(sum(ALPHA,1))) == size(ALPHA,2) );
        end
        
        % Generate M (Sigma) matrix
        M = GAMMA ./ (1 - diag(D));
        % -row sum check
        try
            assert(sum(sum(M,2)) == size(M,1));
        catch
            % Exception: SWE sector 19 row sum equals 0
            sumZeroM = Alpha_text(isnan(sum(M,2)),1:3);
            exception = [yr {'SWE'} {'19'}];
            assert(sum(strcmp(sumZeroM,exception)) == 3);
            
            assert(sum(sum(isnan(M),2)) == n_cs);
            M (isnan(M)) = 0;
            assert(sum(sum(M,2)) == size(M,1) - 1); % 1 less because one of the row totals = 0
        end
               
        % Download sectorwise trade & generate sector weights - SECOND DEFINITION
        SWEIGHTSdata = importdata([sweightsfolder,'WIOD_Sweights_alltimeAverage_secondDef.csv']);
        Sweights_text = SWEIGHTSdata.textdata(2:end,1:1);
        SWEIGHTS = SWEIGHTSdata.data;
        det_sweights = det(SWEIGHTS);
        SWEIGHTS_inv = SWEIGHTS \ eye(size(SWEIGHTS));
        det_sweights(1,2) = det(SWEIGHTS_inv);
        
        % Download supplying-sector-betas
        [num,txt] = xlsread(beta_file,'betas');
        supplying_sector_list = txt(2:end,1);
        BETA_supplying_sector = -1 * num;
        
        % Generate using-sector-betas
        %BETA_using_sector = SWEIGHTS_inv * BETA_supplying_sector; % FIRST DEFINITION
        BETA_using_sector = SWEIGHTS * BETA_supplying_sector;
        
        % Generate B matrix
        B_using_sector = diag(repmat(BETA_using_sector,n_country,1));
        
        % Download total output weights
        TOWEIGHTSdata = importdata([toweightsfolder,'WIOD_TOweights_' yr '.csv']);
        TOweights_text = TOWEIGHTSdata.textdata(2:end,1:3);
        TOWEIGHTS = TOWEIGHTSdata.data;
        
%         fprintf('  --- Imported 17 X 31 matrices \n');
        
        %% Check country - sector alignment across matrices
        
        % Sector list: Gamma versus PPI
        gamma_sector = Gamma_text(:,3);
        ppi_sector   = ppi_text(:,4);
        assert(sum(strcmp(ppi_sector,gamma_sector)) == size(gamma_sector,1));
        
        % Sector list: Gamma versus PPIusd
        ppiusd_sector   = ppiusd_text(:,4);
        assert(sum(strcmp(ppiusd_sector,gamma_sector)) == size(gamma_sector,1));
        
        % Sector list: Sector weights versus Gamma
        assert(sum(strcmp(Sweights_text,gamma_sector(1:n_sector))) == size(Sweights_text,1) );
        
        % Sector list: Sector specific beta versus Gamma
        assert(sum(strcmp(supplying_sector_list,gamma_sector(1:n_sector))) == size(supplying_sector_list,1) );
        B_sector = repmat(supplying_sector_list,n_country,1);
        assert(sum(strcmp(B_sector,gamma_sector)) == size(B_sector,1) );
        
        % Sector list: total output weights versus Gamma
        toweights_sector = TOweights_text(:,3);
        assert(sum(strcmp(toweights_sector,gamma_sector)) == size(gamma_sector,1) );
        
%         fprintf('  --- Sector list: aligned \n');
        
        % Country list: Gamma versus PPI
        gamma_country = Gamma_text(:,2);
        ppi_country = ppi_text(:,5);
        assert(sum(strcmp(ppi_country,gamma_country)) == size(gamma_country,1));
        
        % Country list: Gamma versus PPI
        ppiusd_country = ppiusd_text(:,5);
        assert(sum(strcmp(ppiusd_country,gamma_country)) == size(gamma_country,1));
        
        % Country list: E0 versus Gamma
        E_short_country = E_text(:,4);
        E_country = repelem(E_short_country,n_sector,1);
        assert(sum(strcmp(E_country,gamma_country)) == size(gamma_country,1));
        
        % Country list: total output weights versus Gamma
        toweights_country = TOweights_text(:,2);
        assert(sum(strcmp(toweights_country,gamma_country)) == size(gamma_country,1) );
        
%         fprintf('  --- Country list: aligned for 17 X 31 \n');
        
        %% Generate counterfactual GAMMA matrices
        % Generate autarky GAMMA
        GAMMA_autarky = zeros(size(GAMMA));
        for c = 1:n_country
            for s = 1:n_sector
                line = (c-1)*n_sector+s;
                for ss=1:n_sector
                    %columns of the particular sector in all countries
                    relevant_cols = [ss (1:(n_country-1)).* n_sector + ss];
                    self_column = (c-1)*n_sector + ss;
                    GAMMA_autarky(line,self_column) = sum(GAMMA(line,relevant_cols));
                    %share_self(line,self_column) = GAMMA(line,self_column)/GAMMA_autarky(line,self_column);
                end
            end
        end
        
        
        % Generate smooth world-I GAMMA
        GAMMA_smooth_worldI = zeros(size(GAMMA));
        for c = 1:n_country
            self_lines = (c-1)*n_sector+1:(c-1)*n_sector+n_sector;
            self_columns = self_lines;
            
            foreign_columns = 1:size(GAMMA,2);
            foreign_columns(self_columns) = [];
            
            avg_self = sum(sum(GAMMA(self_lines,self_columns),1),2) / (n_sector*n_sector);
            avg_others = sum(sum(GAMMA(self_lines,foreign_columns),1),2) / (n_sector*n_sector*(n_country-1));
            
            GAMMA_smooth_worldI(self_lines,self_columns) = repmat(avg_self,n_sector,n_sector);
            GAMMA_smooth_worldI(self_lines,foreign_columns) = repmat(avg_others,n_sector,n_sector*(n_country-1));
        end
        
        % Normalize
        GAMMA_sum = sum(GAMMA,2);
        GAMMA_sum_smooth_worldI = sum(GAMMA_smooth_worldI,2);
        ratio = GAMMA_sum ./ GAMMA_sum_smooth_worldI;
        GAMMA_smooth_worldI = GAMMA_smooth_worldI .* repmat(ratio,1,size(GAMMA_smooth_worldI,2));
        
        
        % Generate smooth world-II GAMMA
        GAMMA_smooth_worldII = zeros(size(GAMMA));
        for c = 1:n_country
            c_lines = (c-1)*n_sector+1:(c-1)*n_sector+n_sector;
            for cc=1:n_country
                cc_columns = (cc-1)*n_sector+1:(cc-1)*n_sector+n_sector;
                
                avg_c_cc = sum(sum(GAMMA(c_lines,cc_columns),1),2) / (n_sector*n_sector);
                GAMMA_smooth_worldII(c_lines,cc_columns) = repmat(avg_c_cc,n_sector,n_sector);
            end
        end
        
        % Normalize to have the same share of total input share (so also the same share of value added in total output)
        GAMMA_sum_smooth_worldII = sum(GAMMA_smooth_worldII,2);
        ratio = GAMMA_sum ./ GAMMA_sum_smooth_worldII;
        GAMMA_smooth_worldII = GAMMA_smooth_worldII .* repmat(ratio,1,size(GAMMA_smooth_worldII,2));
        
%         fprintf('  --- Generated counterfactual GAMMAs: autarky, smooth world I-II   \n');
        
        %% --- Generate cost shocks ---

        %% Cost shock: Baseline: beta_cost = beta_fx = 1, beta_comp = 0
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*
        
        % Betas:
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        beta_cost = 1;
        beta_fx   = 1;
        beta_comp = 0;
        
        % BETA matrices
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        BETA_cost = beta_cost * I;
        BETA_fx   = beta_fx   * I;
        BETA_comp = beta_comp * I;
        
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = I;
        LAMBDA_fx   = I;
        % LAMBDAs (_cost & _fx) generated using the formula are very close to the I -Identity matrix,
        % only when rounded to around the 5th decimal point. Therefore used the Identity matrix directly instead.
        % See below for how close it is.
        
        LAMBDA_comp = ZERO;
        
        %{
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = diag(ALPHA' * diag(BETA_cost));
        LAMBDA_fx   = diag(ALPHA' * diag(BETA_fx));
        LAMBDA_comp = diag(ALPHA' * diag(BETA_comp));
        
        % Check if LAMBDA_cost = I & LAMBDA_fx = I & LAMBDA_comp = 0
        assert(round(det(LAMBDA_cost),6) == 1);
        assert(round(det(LAMBDA_fx),6)   == 1);
        assert(round(det(LAMBDA_comp),6) == 0);
        %}
        
        % Compute W
        IminusBd_inv  = (I - BETA_comp) \ I ;
        W_first_term  = (LAMBDA_cost + ALPHA' * IminusBd_inv * BETA_comp * BETA_cost * M) \ I;
        W_second_term = (I - BETA_fx - IminusBd_inv * BETA_comp * BETA_fx);
        W_third_term  = (IminusBd_inv * BETA_comp * BETA_fx * M);
        W = W_first_term * (PPI + ((I - LAMBDA_fx) - ALPHA' * ( W_second_term + W_third_term )) * E_0);
        
        % Check W == PPI
        assert(sum(sum(abs(W -PPI))) / size (W,1) < 1e-08);
        compare_w_ppi = [W PPI];
        
        % Compute cost shock
        COST_baseline = D_inv * ((I - IminusBd_inv * BETA_cost * GAMMA) * W + IminusBd_inv * (BETA_fx * (I - D) - BETA_fx * GAMMA ) * E_0);
        simplified_cost1 = D_inv * ((I - GAMMA) * PPI + (I - D - GAMMA) * E_0);
        
        % Check COST1 = simplified_COST1
        assert(sum(sum(abs(COST_baseline - simplified_cost1))) / size (COST_baseline,1) < 1e-08);
        compare_cost1_cost1simplified = [COST_baseline simplified_cost1];
        
%         fprintf('  --- Generated Cost shocks for baseline  \n');
        
        %% Cost shock: Price complementarities: beta_cost = beta_fx = beta, beta_comp = 1 - beta
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*
        beta_mat = [2/3 1/3]; % 0.33 0.66 out [Andrei, 25.01.2018, use actual 1/3 and 2/3 instead of 0.33 and 0.66]
        
        for b = 1:length(beta_mat)
            beta = beta_mat(b);
            
            % Betas:
            % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
            beta_cost = beta;
            beta_fx   = beta;
            beta_comp = 1-beta;
            
            % BETA matrices
            % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
            BETA_cost = beta_cost * I;
            BETA_fx   = beta_fx   * I;
            BETA_comp = beta_comp * I;
            
            % Generate LAMBDA
            % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
            LAMBDA_cost = beta_cost * I;
            LAMBDA_fx   = beta_fx   * I;
            LAMBDA_comp = beta_comp * I;
            % LAMBDAs (_cost, _fx & _comp) generated using the formula are very close to corresponding betas times I -Identity matrix,
            % only when rounded to around the 5th decimal point. Therefore used the beta times I directly instead.
            % See below for how close it is.
            
            %{
            % Generate LAMBDA
            % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
            LAMBDA_cost = diag(ALPHA' * diag(BETA_cost));
            LAMBDA_fx   = diag(ALPHA' * diag(BETA_fx));
            LAMBDA_comp = diag(ALPHA' * diag(BETA_comp));
            
            % Check if LAMBDA_cost = beta * I & LAMBDA_fx = beta * I & LAMBDA_comp = (1 - beta) * I
            assert(round(det(LAMBDA_cost/beta_cost),6) == 1);
            assert(round(det(LAMBDA_fx/beta_fx),6)     == 1);
            assert(round(det(LAMBDA_comp/beta_comp),6) == 1);
            %}
            
            % Compute W
            IminusBd_inv  = (I - BETA_comp) \ I ;
            W_first_term  = (LAMBDA_cost + ALPHA' * IminusBd_inv * BETA_comp * BETA_cost * M) \ I;
            W_second_term = (I - BETA_fx - IminusBd_inv * BETA_comp * BETA_fx);
            W_third_term  = (IminusBd_inv * BETA_comp * BETA_fx * M);
            W = W_first_term * (PPI + ((I - LAMBDA_fx) - ALPHA' * ( W_second_term + W_third_term )) * E_0);
            
            % Compute cost shock
            COST_pcomp = D_inv * ((I - IminusBd_inv * BETA_cost * GAMMA) * W + IminusBd_inv * (BETA_fx * (I - D) - BETA_fx * GAMMA ) * E_0);
            
            % Cost shocks for different betas
            betaname = replace(num2str(round(beta,3)),'.','');
            eval(['COST_pcomp_' betaname '= COST_pcomp;']);
%             fprintf(['  --- Generated Cost shocks for pcomp, beta = ' betaname '  \n']);
            clear COST_pcomp
        end
        
        %% Cost shock: Price complementarities: sector-specific betas: beta_cost = beta_fx = beta_using_sector, beta_comp = 1 - beta_using_sector
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*
        
        % BETA matrices
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        BETA_cost = B_using_sector;
        BETA_fx   = B_using_sector;
        BETA_comp = I - B_using_sector;
        
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = diag(ALPHA' * diag(BETA_cost));
        LAMBDA_fx   = diag(ALPHA' * diag(BETA_fx));
        LAMBDA_comp = diag(ALPHA' * diag(BETA_comp));
        
        % Check if LAMBDA_cost = LAMBDA_fx & LAMBDA_comp = 0
        assert(sum(sum(LAMBDA_cost - LAMBDA_fx)) == 0);
        assert(round(det(LAMBDA_comp),6) == 0);
        
        % Compute W
        IminusBd_inv  = (I - BETA_comp) \ I ;
        W_first_term  = (LAMBDA_cost + ALPHA' * IminusBd_inv * BETA_comp * BETA_cost * M) \ I;
        W_second_term = (I - BETA_fx - IminusBd_inv * BETA_comp * BETA_fx);
        W_third_term  = (IminusBd_inv * BETA_comp * BETA_fx * M);
        W = W_first_term * (PPI + ((I - LAMBDA_fx) - ALPHA' * ( W_second_term + W_third_term )) * E_0);
        
        % Compute cost shock
        COST_pcomp_ss = D_inv * ((I - IminusBd_inv * BETA_cost * GAMMA) * W + IminusBd_inv * (BETA_fx * (I - D) - BETA_fx * GAMMA ) * E_0);
        
%         fprintf('  --- Generated Cost shocks for pcomp, beta = ss  \n');
        
        %% Cost shock: No FX: E_0 = 0, beta_cost = beta_fx = 1, beta_comp = 0
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*
        
        % Betas:
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        beta_cost = 1;
        beta_fx   = 1;
        beta_comp = 0;
        
        % BETA matrices
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        BETA_cost = beta_cost * I;
        BETA_fx   = beta_fx   * I;
        BETA_comp = beta_comp * I;
        
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = I;
        LAMBDA_fx   = I;
        % LAMBDAs (_cost & _fx) generated using the formula are very close to the I -Identity matrix,
        % only when rounded to around the 5th decimal point. Therefore used the Identity matrix directly instead.
        % See below for how close it is.
        
        LAMBDA_comp = ZERO;
        
        %{
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = diag(ALPHA' * diag(BETA_cost));
        LAMBDA_fx   = diag(ALPHA' * diag(BETA_fx));
        LAMBDA_comp = diag(ALPHA' * diag(BETA_comp));
        
        % Check if LAMBDA_cost = I & LAMBDA_fx = I & LAMBDA_comp = 0
        assert(round(det(LAMBDA_cost),6) == 1);
        assert(round(det(LAMBDA_fx),6)   == 1);
        assert(round(det(LAMBDA_comp),6) == 0);
        %}
        
        % Set elements of E_0 to zero
        E_0_noFX = 0 * E_0;
        
        % Compute W
        IminusBd_inv  = (I - BETA_comp) \ I ;
        W_first_term  = (LAMBDA_cost + ALPHA' * IminusBd_inv * BETA_comp * BETA_cost * M) \ I;
        W_second_term = (I - BETA_fx - IminusBd_inv * BETA_comp * BETA_fx);
        W_third_term  = (IminusBd_inv * BETA_comp * BETA_fx * M);
        W = W_first_term * (PPI + ((I - LAMBDA_fx) - ALPHA' * ( W_second_term + W_third_term )) * E_0_noFX);
        
        % Check W == PPI
        assert(sum(sum(abs(W -PPI))) / size (W,1) < 1e-08);
        compare_w_ppi_noFX = [W PPI];
        
        % Compute cost shock
        COST_noFX = D_inv * ((I - IminusBd_inv * BETA_cost * GAMMA) * W + IminusBd_inv * (BETA_fx * (I - D) - BETA_fx * GAMMA ) * E_0_noFX);
        simplified_cost_noFX = D_inv * ((I - GAMMA) * PPI + (I - D - GAMMA) * E_0_noFX);
        simplified_cost_noFX2 = D_inv * ((I - GAMMA) * PPI );
        
        compare = [COST_noFX simplified_cost_noFX simplified_cost_noFX2];
        
        % Check COST_noE = simplified_COST_noFX
        assert(sum(sum(abs(COST_noFX - simplified_cost_noFX))) / size (COST_noFX,1) < 1e-08);
        compare_costnoFX_costnoFXsimplified = [COST_noFX simplified_cost_noFX];
        
%         fprintf('  --- Generated Cost shocks for noFX   \n');
        
        %% Cost shock: Mechanical pass-through: beta_cost = beta_fx = beta, beta_comp = 0
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*
        beta_mat = [2/3 1/3 ]; % 0.33 0.66 out [Andrei, 25.01.2018, use actual 1/3 and 2/3 instead of 0.33 and 0.66]
        
        for b = 1:length(beta_mat)
            beta = beta_mat(b);
            
            % Betas:
            % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
            beta_cost = beta;
            beta_fx   = beta;
            beta_comp = 0;
            
            % BETA matrices
            % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
            BETA_cost = beta_cost * I;
            BETA_fx   = beta_fx   * I;
            BETA_comp = beta_comp * I;
            
            % Generate LAMBDA
            % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
            LAMBDA_cost = beta_cost * I;
            LAMBDA_fx   = beta_fx   * I;
            LAMBDA_comp   = ZERO;
            % LAMBDAs (_cost&_fx) generated using the formula are very close to corresponding betas times I -Identity matrix,
            % only when rounded to around the 5th decimal point. Therefore used the beta times I directly instead.
            % See below for how close it is.
            
            %{
            % Generate LAMBDA
            % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
            LAMBDA_cost = diag(ALPHA' * diag(BETA_cost));
            LAMBDA_fx   = diag(ALPHA' * diag(BETA_fx));
            LAMBDA_comp = diag(ALPHA' * diag(BETA_comp));
            
            % Check if LAMBDA_cost = beta * I & LAMBDA_fx = beta * I & LAMBDA_comp = 0
            assert(round(det(LAMBDA_cost/beta_cost),6) == 1);
            assert(round(det(LAMBDA_fx/beta_fx),6)     == 1);
            assert(round(det(LAMBDA_comp),6)           == 0);
            %}
            
            % Compute W
            IminusBd_inv  = (I - BETA_comp) \ I ;
            W_first_term  = (LAMBDA_cost + ALPHA' * IminusBd_inv * BETA_comp * BETA_cost * M) \ I;
            W_second_term = (I - BETA_fx - IminusBd_inv * BETA_comp * BETA_fx);
            W_third_term  = (IminusBd_inv * BETA_comp * BETA_fx * M);
            W = W_first_term * (PPI + ((I - LAMBDA_fx) - ALPHA' * ( W_second_term + W_third_term )) * E_0);
            % W = W_first_term * (PPI + ((I - LAMBDA_fx) - ( W_second_term + W_third_term )) * E_0); % ALPHA out
            
            %W_simplified = (BETA_cost \ I) * (PPI + (I - BETA_cost) * (I - ALPHA') * E_0);
            %W_simplified2 = (BETA_cost \ I) * (PPI + (I - BETA_cost) * E_0);
            %www = [W W_simplified W_simplified2];
            
            % Compute cost shock
            COST_mechpt = D_inv * ((I - IminusBd_inv * BETA_cost * GAMMA) * W + IminusBd_inv * (BETA_fx * (I - D) - BETA_fx * GAMMA ) * E_0);
            
            %cost_simplified_mechpt = D_inv * ((I - BETA_cost * GAMMA) * W + BETA_fx * (I - D - GAMMA ) * E_0);
            %cost_simplified_mechpt2 = D_inv * ((I - BETA_cost * GAMMA) * W_simplified2 + BETA_fx * (I - D - GAMMA ) * E_0);
            %www = [www COST_mechpt cost_simplified_mechpt2];
            
            % Cost shocks for different betas
            betaname = replace(num2str(round(beta,3)),'.','');
            eval(['COST_mechpt_' betaname '= COST_mechpt;']);
%             fprintf(['  --- Generated Cost shocks for mechpt, beta = ' betaname '   \n']);
            clear COST_mechpt
        end
        
        %% Cost shock: Mechanical pass-through: sector-specific betas: beta_cost = beta_fx = beta_using_sector, beta_comp = 0
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*
        
        % BETA matrices
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        BETA_cost = B_using_sector;
        BETA_fx   = B_using_sector;
        BETA_comp = ZERO;
        
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = diag(ALPHA' * diag(BETA_cost));
        LAMBDA_fx   = diag(ALPHA' * diag(BETA_fx));
        LAMBDA_comp = diag(ALPHA' * diag(BETA_comp));
        
        % Check if LAMBDA_cost = LAMBDA_fx & LAMBDA_comp = 0
        assert(sum(sum(LAMBDA_cost - LAMBDA_fx)) == 0);
        assert(round(det(LAMBDA_comp),6) == 0);
        
        % Compute W
        IminusBd_inv  = (I - BETA_comp) \ I ;
        W_first_term  = (LAMBDA_cost + ALPHA' * IminusBd_inv * BETA_comp * BETA_cost * M) \ I;
        W_second_term = (I - BETA_fx - IminusBd_inv * BETA_comp * BETA_fx);
        W_third_term  = (IminusBd_inv * BETA_comp * BETA_fx * M);
        W = W_first_term * (PPI + ((I - LAMBDA_fx) - ALPHA' * ( W_second_term + W_third_term )) * E_0);
        
        % Compute cost shock
        COST_mechpt_ss = D_inv * ((I - IminusBd_inv * BETA_cost * GAMMA) * W + IminusBd_inv * (BETA_fx * (I - D) - BETA_fx * GAMMA ) * E_0);
        
%         fprintf('  --- Generated Cost shocks for mechpt, beta = ss  \n');

        %% Counterfactual PPIs: Input linkages: Balanced I & II, domestic linkages: beta_cost = beta_fx = 1, beta_comp = 0
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*

        % BALANCED I
        % Use baseline cost shock
        cost_balI = COST_baseline;
        
        % Compute counterfactual PPI
        % - Replace the GAMMA with GAMMA_smooth_worldII % note the twist in the name 
        IminusGamma_inv  = (I - GAMMA_smooth_worldII) \ I ;
        IminusDminusGamma = I - D - GAMMA_smooth_worldII;
        PPI_balI = IminusGamma_inv * (D* cost_balI - IminusDminusGamma * E_0);
        
%         fprintf('  --- Generated Counterfactual PPI -balanced I-  \n');
        
        % BALANCED II
        % Use baseline cost shock
        cost_balII = COST_baseline;
        
        % Compute counterfactual PPI
        % - Replace the GAMMA with GAMMA_smooth_worldI % note the twist in the name 
        IminusGamma_inv  = (I - GAMMA_smooth_worldI) \ I ;
        IminusDminusGamma = I - D - GAMMA_smooth_worldI;
        PPI_balII = IminusGamma_inv * (D * cost_balII - IminusDminusGamma * E_0);
        
%         fprintf('  --- Generated Counterfactual PPI -balanced II-  \n');
        
        % DOMESTIC LINKAGES
        % Use cost shock no FX
        cost_domLin = D * COST_baseline; % ? no divide: no D^-1
        
        % Compute counterfactual PPI
        % - Replace the GAMMA with GAMMA_autarky
        IminusGamma_inv  = (I - GAMMA_autarky) \ I ;
        PPI_domLin = IminusGamma_inv * cost_domLin; % ? no divide: no D^-1
        
%         fprintf('  --- Generated Counterfactual PPI -domestic linkages-  \n');
        
        % Check PPI = counterfactual PPIs
        compare_ppi_ppicounterfactual = [PPI PPI_balI PPI_balII PPI_domLin] ;
%         plot(compare_ppi_ppicounterfactual,'DisplayName','compare_ppi_ppicounterfactual')

%         corr(compare_ppi_ppicounterfactual)
        
%         fprintf('  --- Correlation of PPI and counterfactual PPIs- linkages \n');
        
        %% Counterfactual PPI: No oil
        selector_noOil = 1 - strcmp(Gamma_text(:,3),'23'); % energy sector
        Shock_noOil = selector_noOil .* COST_baseline;
        IminusGamma_inv = (I - GAMMA) \ I ;        
        PPI_noOil = IminusGamma_inv * D * Shock_noOil;
        
%         fprintf('  --- Generated Counterfactual PPI -no oil-  \n');
        
        % Check PPI = counterfactual PPIs
        compare_ppi_ppinoOil = [PPI PPI_noOil] ;
%         plot(compare_ppi_ppinoOil,'DisplayName','compare_ppi_ppinoOil')

%         corr(compare_ppi_ppinoOil)
        
%         fprintf('  --- Correlation of PPI and counterfactual PPI- no oil \n');
        
        %% Cost shock: Baseline with USD denominated PPI: beta_cost = beta_fx = 1, beta_comp = 0
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*
        
        % Betas:
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        beta_cost = 1;
        beta_fx   = 1;
        beta_comp = 0;
        
        % BETA matrices
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        BETA_cost = beta_cost * I;
        BETA_fx   = beta_fx   * I;
        BETA_comp = beta_comp * I;
        
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = I;
        LAMBDA_fx   = I;
        % LAMBDAs (_cost & _fx) generated using the formula are very close to the I -Identity matrix,
        % only when rounded to around the 5th decimal point. Therefore used the Identity matrix directly instead.
        % See below for how close it is.
        
        LAMBDA_comp = ZERO;
        
        %{
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = diag(ALPHA' * diag(BETA_cost));
        LAMBDA_fx   = diag(ALPHA' * diag(BETA_fx));
        LAMBDA_comp = diag(ALPHA' * diag(BETA_comp));
        
        % Check if LAMBDA_cost = I & LAMBDA_fx = I & LAMBDA_comp = 0
        assert(round(det(LAMBDA_cost),6) == 1);
        assert(round(det(LAMBDA_fx),6)   == 1);
        assert(round(det(LAMBDA_comp),6) == 0);
        %}
        
        % Set elements of E_0 to zero
        E_0_noFX = 0 * E_0;
        
        % Compute W
        IminusBd_inv  = (I - BETA_comp) \ I ;
        W_first_term  = (LAMBDA_cost + ALPHA' * IminusBd_inv * BETA_comp * BETA_cost * M) \ I;
        W_second_term = (I - BETA_fx - IminusBd_inv * BETA_comp * BETA_fx);
        W_third_term  = (IminusBd_inv * BETA_comp * BETA_fx * M);
        W = W_first_term * (PPIusd + ((I - LAMBDA_fx) - ALPHA' * ( W_second_term + W_third_term )) * E_0_noFX);
        
        % Check W == PPI
        assert(sum(sum(abs(W -PPIusd))) / size (W,1) < 1e-08);
        compare_w_ppi = [W PPIusd];
        
        % Compute cost shock
        COST_ppiusd = D_inv * ((I - IminusBd_inv * BETA_cost * GAMMA) * W + IminusBd_inv * (BETA_fx * (I - D) - BETA_fx * GAMMA ) * E_0_noFX);
        simplified_costppiusd = D_inv * ((I - GAMMA) * PPIusd + (I - D - GAMMA) * E_0_noFX);
        
        % Check COST1 = simplified_COST1
        assert(sum(sum(abs(COST_ppiusd - simplified_costppiusd))) / size (COST_ppiusd,1) < 1e-08);
        compare_cost1_cost1simplified = [COST_ppiusd simplified_costppiusd];
        
%         fprintf('  --- Generated Cost shocks for baseline with USD denominated PPI   \n');
        
        %% --- Prepare output ---
        %% Gather PPI & cost shocks in a table
        headers_textdata = ppidata.textdata(1,1:5);
        textdata = ppidata.textdata(2:end,1:5);
        
        % Get the list of all generated cost shocks, and add PPI
        if y == 1995 && m == 1
            numdatanames = who('-regexp','PPI*');
            numdatanames = [numdatanames; who('-regexp','COST*')];
        end
        
        numdata = []; headers_numdata = ''; % initialize
        for n = 1:length(numdatanames)
            numdata = [numdata eval(numdatanames{n})];
            headers_numdata = [headers_numdata cellstr(numdatanames{n})];
        end
        
        %headers_numdata = ['PPI' 'COST1' 'COST1_simplified' 'COST2_033' 'COST2_067' 'COST3_033' 'COST3_067' 'COST4' 'COST5' ];
        %numdata = [PPI COST1 COST1_simplified COST2_033 COST2_067 COST3_033 COST3_067 COST4 COST5 ];
        
        TABLE = [];
        for i = 1: length(headers_textdata)
            clear T;
            T = table(textdata(:,i),'variableNames',cellstr(headers_textdata{i}));
            TABLE = [TABLE T];
        end
        
        for i = 1: length(headers_numdata)
            clear T;
            T = table(numdata(:,i),'variableNames',cellstr(headers_numdata{i}));
            TABLE = [TABLE T];
        end
        
        %% Generate aggregated PPI & cost shocks at the country level
        % - using 2002 total output weights
        
        % Download total output weights 2002
        TOWEIGHTSdata2002 = importdata([toweightsfolder,'WIOD_TOweights_2002.csv']);
        TOweights2002_text = TOWEIGHTSdata2002.textdata(2:end,1:3);
        TOWEIGHTS2002 = TOWEIGHTSdata2002.data;
        TOWEIGHTS2002_tomultiply = reshape(TOWEIGHTS2002,[n_sector,n_country]);
        
        % Reshape for multiplying by weighting matrix
        for n = 1:length(numdatanames)
            eval([ numdatanames{n} '_forweighing = transpose(reshape(' numdatanames{n} ',[n_sector,n_country]));']);
        end
        
        % Multiply and take the diagonal: to weigh
        for n = 1:length(numdatanames)
            eval([ numdatanames{n} '_ag = diag(' numdatanames{n} '_forweighing * TOWEIGHTS2002_tomultiply);']);
        end
        %% Gather aggregated PPI & cost shocks in a table
        
        % Get the list of generated aggregated cost shocks, and add PPI
        numdata_aggregate = []; headers_numdata_aggregate = ''; % initialize
        for n = 1:length(numdatanames)
            numdata_aggregate = [numdata_aggregate eval([numdatanames{n} '_ag'])];
            headers_numdata_aggregate = [headers_numdata_aggregate cellstr([numdatanames{n} '_ag'])];
        end
        
        %headers_numdata_aggregate = list2Cell('PPI_aggregate,COST1_aggregate');
        %numdata_aggregate = [PPI_aggregate COST1_aggregate];
        
        % Prepare text data
        selectC = [1:n_sector:n_cs];
        country_text = Gamma_text(selectC,2);
        year_text = cellstr(repmat(yr,size(country_text)));
        month_text = cellstr(repmat(mo,size(country_text)));
        
        % Generate table
        Tcountry = table(country_text,'variableNames',{'country'});
        Tyear    = table(year_text,'variableNames',{'year'});
        Tmonth   = table(month_text,'variableNames',{'month'});
        AGGTABLE = [Tcountry Tyear Tmonth];
        
        for i = 1: length(headers_numdata_aggregate)
            clear T;
            T = table(numdata_aggregate(:,i),'variableNames',cellstr(headers_numdata_aggregate{i}));
            AGGTABLE = [AGGTABLE T];
        end
        
        %% Stack tables for each period in a large table with all periods
        if y == 1995 && m == 1
            TABLE_allP = []; AGGTABLE_allP = []; % initialize
        end
        TABLE_allP = [TABLE_allP; TABLE]; % cs breakdown
        AGGTABLE_allP = [AGGTABLE_allP; AGGTABLE]; % c breakdown
        
        %% -- (End the loop for the month)
    end
    
    %% Stack BETA_using_sector into one long vector
    if y == 1995 && m == 12
        BETA_using_sector_allP = []; % initialize
    end
    BETA_using_sector_allP = [BETA_using_sector_allP; BETA_using_sector];
    
    %% -- (End the loop for the year)
end


%% Export to excel/csv

% Beta_u, Beta_s and S
output_filename = 'betas_table.xlsx';
num_output_beta = [BETA_supplying_sector BETA_using_sector SWEIGHTS];
text_output_beta = supplying_sector_list;
headers_beta = ['Sector' 'Beta_supplying_sector' 'Beta_using_sector' supplying_sector_list'];
output_beta = [text_output_beta num2cell(num_output_beta)];
output_beta = [headers_beta; output_beta];
TABLE_beta = table(output_beta);
writetable(TABLE_beta,fullfile(outputfolder,output_filename),'WriteVariableNames',false);
fprintf(['  --- Betas (sector specific) exported to excel: ' output_filename '  \n']);

% Cost shock for each country-sector
output_filename = ['costShocks_allPeriods_table.xlsx'];
writetable(TABLE_allP,fullfile(outputfolder,output_filename),'WriteRowNames',true);
fprintf(['  --- Cost shocks exported to excel: ' output_filename '  \n']);

% Aggregated cost shock for each country
output_filename = ['aggregatedCostShocks_allPeriods_table.xlsx']; % The file used in factor analysis
writetable(AGGTABLE_allP,fullfile(outputfolder,output_filename),'WriteRowNames',true);
fprintf(['  --- Aggregated cost shocks exported to excel: ' output_filename '  \n']);

toc
pause(3)
exit
end


